#include <iostream>
#include "Spline.h"
#include "iofile.h"
#include <cmath>

class X :public Function{
public:
    double operator() (const double &st) const{
        return sqrt(3)*sin(st);
    }

    double diff(const double &st) const{
        return sqrt(3)*cos(st);
    }

    double diff2(const double &st) const{
        return -sqrt(3)*sin(st);
    }
};

class Y :public Function{
public:
    double operator() (const double &st) const{
        return 2.0/3.0*(  sqrt(3)*cos(st) + sqrt(  sqrt(3)*abs(sin(st)) ));
    }
    double diff(const double &st) const{
        double tmp = ((*this)(st+tao)-(*this)(st-tao)) / (2*tao);
        if(tmp>10) tmp=10;
        else if(tmp<-10) tmp=-10;

        return tmp;
    }

    double diff2(const double &st) const{
        double tmp =((*this)(st+2*tao)+(*this)(st-2*tao)-2*(*this)(st)) / (4*tao*tao);
        if(tmp>10) tmp=10;
        else if(tmp<-10) tmp=-10;
        return tmp;
    }
};

class Heart{
public:
    X xs;
    Y ys;
    double L = 2.0*pi;
    double curve_S ;
    int M;
    BCT bct;
    vector <double> bct_x;
    vector <double> bct_y;
    std::vector<double> s;
    std::vector<double> knot_xs;
    std::vector<double> knot_ys;

public:
    //生产knots
    Heart(int knots_num ,BCT bct):bct(bct){
        //设置边界条件
        bct_x.clear();
        bct_y.clear();
        switch (bct){
        case complete:
            bct_x.push_back(xs.diff(0));
            bct_x.push_back(xs.diff(L));
            bct_y.push_back(ys.diff(0));
            bct_y.push_back(ys.diff(0));
        break;

        case cubic:
            bct_x.push_back(xs.diff2(0));
            bct_x.push_back(xs.diff2(L));

            bct_y.push_back(ys.diff2(0));
            bct_y.push_back(ys.diff2(0));
        break;

        case natural:
            bct_x.push_back(0);
            bct_x.push_back(0);

            bct_y.push_back(0);
            bct_y.push_back(0);
        break;
        default:
            break;
        }

        //设置knots
        M = int(knots_num/2.0+1);
        int m = int(M/6.0); 
        double eps_singu_pt =pi/10.0;
        if(M>15) eps_singu_pt = pi/20.0;

        for(int i=0;i<=m;i++){
            double tmp = eps_singu_pt*i/double(m);
            s.push_back(tmp);
            knot_xs.push_back(xs(tmp));
            knot_ys.push_back(ys(tmp));
        }

        for(int i=1;i<=M-2*m;i++){
            double tmp = eps_singu_pt+
                        (L/2.0-2.0*eps_singu_pt)*i/double(M-2*m);
            s.push_back(tmp);
            knot_xs.push_back(xs(tmp));
            knot_ys.push_back(ys(tmp));
        }

        for(int i=1;i<=m;i++){
            double tmp = (L/2.0-eps_singu_pt)+eps_singu_pt*i/double(m);
            s.push_back(tmp);
            knot_xs.push_back(xs(tmp));
            knot_ys.push_back(ys(tmp));
        }
        int tmp_size = s.size();

        for(int i=1;i<=tmp_size-1;i++){
            double tmp =2*pi-s[tmp_size-1-i];
            s.push_back(tmp);
            knot_xs.push_back(-knot_xs[tmp_size-1-i]);
            knot_ys.push_back(knot_ys[tmp_size-1-i]);
        }
    }


    //拟合曲线
    void Inter_fitting(){
        vector<vector <double >> knots;
        knots.push_back(knot_xs);
        knots.push_back(knot_ys);

        vector<vector <double >> bcts;
        bcts.push_back(bct_x);
        bcts.push_back(bct_y);
        Curve<2> c(knots,bcts,bct);
        vector<double> Bx;
        vector<double> By;

        //拟合点
        int fitnum = 1000;
        curve_S = c.Curve_S;
        for(int j=0; j<= fitnum;j++){
            Bx.push_back(c.ppspline_fit[0](j * curve_S/float(fitnum)));
            By.push_back(c.ppspline_fit[1](j * curve_S/float(fitnum)));
        }

        ClearFile("E.txt");
        VectorToFile(Bx , "E.txt");
        VectorToFile(By , "E.txt");
        VectorToFile(knot_xs , "E.txt");
        VectorToFile(knot_ys , "E.txt");
    }

};